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Abstract 

The aim of this note is to present a numerical method to solve the Stokes problem in a bounded 
domain with a Dirac source term, which preserves optimality for any approximation order by the 
finite element method. It is based on the knowledge of a fundamental solution of the associated 
operator over the whole space. This method is motivated by the modeling of the movement of active 
thin structures in a viscous fluid. 
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Resume 

Une methode numerique pour la resolution du probleme de Stokes avec une force 
ponctuelle en terme source. Le but de cette note est de presenter une methode numerique pour 
la resolution du probleme de Stokes avec une force ponctuelle en terme source, qui assure l’optimalite 
de l’erreur d’approximation elements finis. Elle s’appuie sur la connaissance explicite d’une solution 
fondamentale de l’operateur lineaire associe. Cette methode est motivee par la nrodelisation du 
mouvenrent de structures fines actives dans un fluide visqueux. 

Mots cles : estimations d’erreur, methode elements finis, Stokeslet, structures fines. 


Version francaise abregee 

L’etude du nrouvement de structures fines actives dans un fluide visqueux, tels que les flagelles 
permettant la nage de bacteries ou les cils impliques dans le transport mucociliaire, conduit a 
considerer le probleme de Stokes avec un second membre singulier. Dans l’asymptotique d’un cil 
dont le diametre tend vers 0 et la vitesse vers l’infini, le terme source est en fait une distribution 
lineique de forces. Dans le but de pouvoir faire des calculs, puisque integrer numeriquement le long 
d’une courbe quelconque est difficile, nous approchons la distribution lineique de forces dr par une 
sonmre de forces ponctuelles ^ Cjdj. Une preuve basee sur celle du theoreme des sommes de Riemann 
permet de montrer qu’il y a convergence, au sens faible dans iW 3 ' 2_s , pour tout s > 0, de ^ c i$i vers 
dr lorsque le nombre N de masses de Dirac tend vers l’infini. On peut aussi preciser la convergence 
dans des espaces plus faibles, voir (1). La convergence des solutions associees se deduit de l’inegalite 
(2), tiree de [1]. On est alors rarnene a l’etude du probleme de Stokes avec une force ponctuelle en 
terme source. 

Lorsqu’on considere un probleme elliptique avec une masse de Dirac en second membre, en 
dimension d ^ 2, ce second membre n’etant pas dans H , le probleme sort du cadre variationnel 
standard base sur l’espace de Sobolev H 1 . Si la methode des elements finis peut etre definie au 
niveau discret, les resultats de convergence classiques ne sont a priori plus valables. Dans le cas 
du probleme de Poisson, qui peut etre vu comme une version scalaire et simplifiee du probleme de 
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Stokes, Scott a demontre dans [ ] que la methode elements finis Pi converge en norme L 2 a l’ordre 1 
en 2d et 1/2 en 3d. Des estimations similaires ont ete obtenues dans [3] avec une methode de Galerkin 
discrete. De plus, Apel et ses co-auteurs ont montre dans [: ] qu’en raffinant le maillage autour de 
la singularite, on retrouvait l’ordre de convergence classique. La methode presentee, basee sur la 
connaissance explicite d’une solution fondamentale de l’operateur lineaire associe, fait partie d’une 
classe de methodes dites de soustraction, introduites en electroencephalographie [ ]. Elle perrnet de 
retrouver les ordres de convergence classiques sans raffinement de maillage. 

Pour fixer les idees, nous allons nous interesser au probleme de Stokes avec des conditions aux 
limites de type Dirichlet homogenes, voir le probleme (4). La particularity de ce probleme reside en 
la singularite du second membre : un Dirac de force applique en un point xq du domaine fi. Pour cet 
operateur, on commit une solution fondamentale definie en domaine infini, appelee Stokeslet, que 
Ton note (u^,p^), voir (5). On obtient la solution (u ,p) du probleme (4) en ajoutant a (u^,^) un 
relevement regulier prenant ainsi en compte les conditions aux bords. La singularite de la solution 
(u ,p) est contenue dans la solution fondamentale (us,ps), et elle est localisee au point xo- Le 
principe de la methode qui suit, est de capturer cette singularite pour se ramener a la resolution 
d’un probleme auxiliaire regulier. 

On commence par definir une fonction plateau x, reguliere, valant 1 sur un voisinage de xo et 
0 loin de ce point, voir Definition 1. On note ensuite uo = x u <5 et po = XP5-, et g et h les fonctions 
definies en ( 6 ). D’apres ces definitions, on remarque que les supports de g et h sont contenus dans 
une couronne centree en xo, voir Figure 1. De plus, les fonctions 11,5 et ps etant analytiques en dehors 
de xo, la regularity des fonctions g et h depend directement de celle de la fonction x- Finalement, 
pour obtenir la solution (u,p) de (4), il suffit de corriger les termes d’erreur g et h introduits en 
( 6 ) en resolvant le probleme elliptique regulier (7), dont on note v la solution. En effet, la fonction 
u := uo + v est la solution du probleme (4). 

Cette methode perrnet de passer de la resolution d’un probleme singulier a celle d’un probleme 
auxiliaire regulier. Alors que le premier converge a un ordre faible [ ], le second converge a l’ordre 
optimal, quel que soit l’ordre des elements utilises. En notant u/j := uo + v/>, ou est la solution 
numerique du probleme (7) obtenue par une methode elements finis, on deduit de ( 8 ) que l’erreur 
commise sur u est la meme que celle commise sur v, et on montre ainsi que la vitesse de convergence 
est optimale. 

Par exemple, si on utilise une methode elements finis P\, k = 0 suffit. On definit alors x comme 
en (9), et on explicite g et h, valant respectivement (10) et (11) en dimension 2, et (12) et (13) 
en dimension 3. Apres resolution numerique du probleme (7), on obtient finalement une solution 
approchee dont l’erreur |u — u ^|| L 2 est en 0 (/r 2 ), quelle que soit la dimension, contre une erreur, 
avec une methode directe, en 0(h ) en dimension 2 et en 0(\/h) en dimension 3. 

Cette methode, presentee dans le cas du probleme de Stokes, peut se generaliser a d’autres 
problemes elliptiques lineaires, comme le probleme de Poisson avec une masse de Dirac en second 
membre. Les conditions aux limites de type Dirichlet homogenes peuvent aussi etre remplacees par 
des conditions de type Dirichlet non homogenes, Neumann ou Robin. Enfin, la linearite, qui joue un 
role essentiel, perrnet en outre de resoudre le cas ou le second membre est la somme d’un nornbre fini 
de forces ponctuelles et d’une fonction lisse, tout en ne resolvant qu’un seul probleme numerique. 

1 Introduction. 

In order to model active thin structures in a viscous fluid, such as flagella connected to bacteria 
or cilia involved in the mucociliary transport, we have studied the Stokes problem with a singular 
right-hand side. In the asymptotic of a zero diameter cilia with an infinite velocity, the source term 
is a lineic distribution of forces, which, in order to ease computations, will be approximated by 
a sum of punctual forces. After having justified this approximation, we will present a numerical 
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method to solve the Stokes problem with a punctual force in source term, and illustrate the results 
by numerical simulations. 


2 Approximation of the lineic distribution of forces by a sum of 
punctual forces. 


Since calculating an integral on any curve is numerically very difficult, the source term, noted 5r, 
the lineic distribution of forces on a curve F, is approached by a sum of N punctual forces 
uniformly distributed along T. The theorem of Riemann sums ensures that abi weakly converges 
to hr in iF~ 3 / 2 ~ s , for all s > 0. Working in weaker spaces, it is possible to adapt the proof of 
theorem of Riemann sums and specify the convergence : 


N 

A - 2 cA 

i=l 


H~ 2 —* 


c 

Vn 


and 


N 

A - cA 

i =1 


H~ 5 / 2 —^ 



(1) 


Moreover, using a result proved by Lions and Magenes in [4], which can be written in this case 


^11 # 2 — r ^ C*|| f || JJ—r , Vr ^ 0, 


( 2 ) 


where u is the solution of a regular elliptic problem with a source term f e H ~ r , we can conclude 
that the solution ujy of the Stokes problem with ]F c,5j right-hand side converges to the solution ur 
of the Stokes problem with hr source term, when N goes to infinity. Actually, we have 


«r _ mjv 



and 


C 

Ur - UN II—1./2—s < 


( 3 ) 


Finally, the solution of the Stokes problem with a lineic distribution of forces is approached by the 
solution of Stokes problem with a finite sum of punctual forces in source term. By linearity and 
without loss of generality, in the following we will deal with a single punctual force. 


3 Numerical method to solve the Stokes problem with a Dirac 
source term. 

In dimension d ^ 2, the 5-distribution is not continuous on H 1 , and so the solution of an elliptic 
problem with Dirac source term is not regular. Consequently, classical results for the convergence 
of the finite element method are not valid. In the case of the Poisson problem, which can be seen as 
the scalar version of the Stokes problem, Scott has shown in [ ] that the Pi-finite element method 
converges for L 2 -norm at the order 1 in dimension 2 and at the order 1/2 in dimension 3. Similar 
estimates have been obtained in [3] with a discrete Galerkin method. Moreover, it has been shown 
by Apel and his co-authors [2] that using graded meshes, it is possible to get numerically the classical 
order of convergence. The aim of this section is to present a numerical method which preserves 
optimality for any approximation order, without using mesh grading. It is based on the knowledge 
of a fundamental solution of the considered linear elliptic problem. This approach fits on the class 
of subtraction methods , introduced in [ ] in the context of electroencephalography. 


3.1 Principle of the method. 

Let us consider the following problem, defined on a bounded open domain P c M d , 


—/jAu + Vp = 
div u = 


Sxo F 

0 

0 


u 
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in P, 
in P, 
on 5P, 


( 4 ) 












where xq is fixed in fl and F is a vector of M. d . Let us note that a fundamental solution of problem 


(4) is known in dimensions 2 and 3 : 




• d = 2, 

M*)-4L(-MM)i 2 + 0)f 

and 

Ps(x) = 

J_xF 

27r |x 2 

• d = 3, 

"iW - iL (Sr + wO F - 

and 

Ps(x) = 

J_xF 

47 r ]xp 


where 1^ is the identity matrix. The fundamental solution (u g,Ps) does not satisfy the boundary 
conditions, and so it is not the solution of problem (4). But this solution can be retrieved by adding 
a regular lifting, therefore the whole information on the singularity of the solution (u ,p) is contained 
in the fundamental solution (u s,ps) and is located at xq. In order to extract this singularity, let us 
fix 0 < a < b < d(x o, dfl) and define X by Definition 1. 


Definition 1. Assume that X is a bump 
function satistying for some k ^ 0, 

• X eH 2 +k (R d ), 

X l B(x 0 ,a) 

* Vx„,6r = 0 - 

Figure 1: Domain D. 

Figure 1: Domaine D. 

Then, with uo := and po := X PS, we define g and h as 

g = —/iAuo + Vpo — <5x 0 F and h = div uo- ( 6 ) 

By the definitions of u, 5 , p ,5 and X , supp( g) c 1Z b a (x 0 ) and supp(h) c; lZ b a (x 0 ), where TZ b (x 0 ) is the 
ring centered around xq, of internal radius a and external radius 6, see Figure 1. Since U 5 and ps 
are analytic on D\{xo}, the regularity of functions g and h directly depends on the regularity of 
function X , namely g e H k (£l) and h e H k+1 (Q). Finally, it only remains to correct the terms uo 
and po by solving the regular elliptic problem 

-gAv + Vg = —g in D, 

div v = —h in D, (7) 

v = 0 on <?D, 

and the solution of problem (4) is given by (u ,p) := (uo + v,po + o) = (x u <5 + v > XPS + q), where uo 
and po are explicitly known functions and (y,q) is the solution of problem (7). Noting (yhiQh.) the 
numerical solution of problem (7) and defining := + uo and ph = qu + Po, we have, 

ll u - Ufc||tf S (fi) = IIv - v ft ||^ S (Q), forO^s^fc + l, 
ip-Ph\\H°(n) = \\q-qh\\H*(n), iorO^s^k. 

Actually, this method allows us to switch from the numerical computation of the solution of a sin¬ 
gular problem with Dirac source term (with a poor convergence rate) to the numerical computation 
of the solution of a regular problem with an optimal convergence rate, at any required precision in 
terms of regularity. 
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3.2 Practical aspects. 

For the sake of simplicity, the location of the Dirac source term will be the origin. First, we need to 
choose a suitable function y. Actually, to take advantage of using P^-hnite elements, £ ^ 1, x has 
to be iP +1 (M d ), in order to ensure that g e I7^ _1 (f2) and h e H e (Q), and finally to get an optimal 
order of convergence. For instance, for £ = 1, let us define y, as a radial function, by: 


x(0 


r 1 

2r 3 — 3 (a + b)r 2 + 6 abr + b 2 (b — 3a) 
(b — o ) 3 

l o 


for r e [0, a], 

for r e [a, b ], 
for r > b, 


(9) 


where the function r is defined on M. d by 


r(xj = ||x|| 2 . 

The function g and h can be explicited. According to this definition of y, g and h vanish outside 
the ring a < ||x ||2 < b. For a < ||x|| 2 < b , the expressions of g and h depend on the dimension, 

• for d = 2, 


g( x ) = 


27r(6 — a) 3 r 


(3r 2 — 2 (a + b)r + ab ) In r + 2 r 2 — 2 (a + b)r + 2 ab^j I 2 + (ab — r 2 ) 


x 4 x 


, . . 3(1 —In r)(r 2 — (a + b)r + ab) ^ 

' !(X) " - 2 Mb-a)^r - X ' F ' 


( 10 ) 

( 11 ) 


for d = 3, 


g(x) = 


/i(x) = 


47 r(6 — a) 3 r 2 

3(r 2 — (a + b)r + ab) 
271 — a) 3 r 2 


((a + b)r — 2 r 2 )l 3 + (2 ab — (a + b)r) 


x*x 


x F. 


( 12 ) 

(13) 


4 Numerical illustrations. 

In this section, we illustrate our theorical results by a numerical example. We define P as the 
unit square and xq = (0.5,0.5). The following table presents the L 2 -error for a direct method (dir. 
meth.) and a subtraction method (sub. meth.) respectively, for a characteristic mesh size h, and 
the estimated order of convergence (e.o.c.). Figure 2 illustrates the section {y = 0.5} of the error 
|u — u/i| in the both cases. Numerical simulations evidence the fact that solving the auxiliary prob¬ 
lem associated to the subtraction procedure of the singularity is more efficient than solving directly 
the problem with the Dirac source term. 


h 

2 ~ a 

2 -4 

2 b 

2 -b 

2 = ‘ 

e.o.c. 

Dir. meth. 

1.02 x 10~ 2 

4.87 x 10~ 3 

2.36 x 10~ 3 

1.21 x 10" 3 

5.89 x 10~ 4 

1.02 

Sub. meth. 

4.12 x 10~ 3 

1.33 x 10" 3 

2.92 x 10" 4 

6.86 x 10" 5 

2.71 x 10" 5 

1.88 
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Figure 2: Section {y = 0.5} of the error |u — u^| for the direct method and the subtraction method 
with h = 0.125. 

Figure 2: Coupe {y = 0.5} de l’erreur |u —u^| pour la methode directe et la methode de soustraction 
avec h = 0.125. 

5 Conclusion. 

To model active thin structures in a viscous fluid, such as flagella connected to bacteria or cilia 
involved in the mucociliary transport, we have studied Stokes problem with a singular right-hand 
side: a punctual force. However, when this problem is solved numerically, the singularity causes 
a poor convergence of the approximate solution to the exact solution. The method presented in 
this note preserves optimality for any approximation order, without using mesh grading. If the 
examples are treated with homogeneous Dirichlet conditions, the same method is still valid in the 
case of non homogenous Dirichlet or any affine boundary conditions (Neumann, Robin...), up to 
suitable adaptations. Similarly, the method can be generalized to the problem with a sum of a finite 
number of Dirac masses and a smooth term right-hand side. 
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